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Abstract 

In recent years more and more long-term broadband data sets are collected 
in geosciences. Therefore there is an urgent need of algorithms which semi- 
automatically analyse and decompose these data into separate periods which 
are associated with different processes. Often Fourier and Wavelet Trans- 
form is used to decompose the data into short and long period effects but 
these fail often because of their simplicity. In this paper we investigate the 
novel approaches Empircial Mode Decomposition and Sparse Decomposition 
for long-term sea floor pressure data analysis and compare them with the 
classical ones. Our results indicate that none of the methods fulfils all the 
requirements but Sparse Decomposition performed best except for computing 
efficiency. 

Keywords: Time Series, Sea Floor Pressure Data, Fourier Transform, 
Wavelet Transform, Empirical Mode Decomposition, Sparse Decomposition 



1. Introduction 
1.1. Motivation 

Long-term monitoring of earth processes becomes more and more essen- 
tial and indispensable as only long records will reveal small changes which 
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may eventually have large societal impact in the near future. These pro- 
cesses encompasses geohazards like monitoring of potential landslides and 
volcanic activity but also processes related to earthquake precursors and cli- 
matic signals. In addition real-time observations of tectonic processes with 
high-precision GPS measurements and satellite based observations of the 
magnetic and gravity field of the earth give new insights into geological pro- 
cesses at local, regional and plate-size scales. As these signals are mostly 
small in amplitude and often buried in large-amplitude signals caused by 
other sources or only visible after major corrections to the field observations 
have been made it is a major challenge to separate the wanted from the 
unwanted signal components. 

Broadband data acquisition is now very common as data storage ca- 
pacity is no major obstacle today and low power data acquisition with au- 
tonomous monitoring systems over considerably long time periods (i.e. years) 
in remote places like the deep ocean floor become more and more feasible. 
In addition large observational networks - either regionally only temporary 
like EarthScope (http://www.earthscope.org/) - or permanent like Nep- 
tune Canada (http : / /www . neptunecanada . ca/) or national or international 
geodetic GPS observation networks deliver continuously time series whose 
information content is not completely taken advantage of. All time series 
contain signal periods from short (several seconds) to long (several month 
to years); origin and cause of some signal components is sometimes quite 
well known (deterministic) as e.g. tides from various origins but in a lot of 
cases especially long-term changes are not well understood. In addition in- 
strumental behaviour like drift and effects of environmental changes severely 
alter the signal. 

In earth science seismology has faced this problem for a long time which 
became even more pressing over the last decade during which the number 
of seismological stations increased accompanied by a tremendous increase of 
digital storage capacity and a rapid decrease in storage cost. In seismology 
robust methods were established to extract earthquakes from long data sets 
automatically and also detect phases and locate earthquakes. This highly au- 
tomated processing can be achieved because the characteristics of the signals 
are well-known and years of research have been spend until these algorithms 
were developed (Bormann, 2002). 

Traditional tools for data analysis are mostly based on spectral analysis 
methods which are well established and available with a lot of software tools. 
However they tend to introduce unwanted effects especially at short periods 
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due to filters and finite length of time series. Also their requirement of a 
stationary time series is often not fulfilled in the case of long-term records. 
Therefore new data exploration tools are needed which enable us to separate 
the different signal components. These tools should be capable of analysing 
long time series efficiently and help to identify different signal contributions 
originating from different processes. 

In this paper we present the results of time series analysis using four dif- 
ferent methods. The times series analysed are sea floor pressure records from 
two different locations (details see section 1.2). The motivation for this in- 
vestigation came from the fact that the sea floor pressure signal is dominated 
by a large deterministic tidal pressure signal but due to the high pressure 
resolution and a technically possible sampling interval of down to 2 seconds, 
small and short transient signals can be observed. Superimposed on these 
tides are also long-period transient changes, associated with tectonic and/or 
oceanographic phenomena but also short period signals originating from local 
or teleseismic events. The problem was and still is to separate the determin- 
istic from the small-amplitude non-deterministic signals. Traditional tidal 
analysis methods as described e.g. in Pawlowicz et al. (2002) and available 
as Matlab® toolbox (http://www.eos.ubc.ca/~rich/) do not always fulfil 
the requirements. 

We tested four methods in total: the traditional (1) Harmonic Decomposi- 
tion, (2) Wavelet Decomposition, (3) Empirical Mode Decomposition (EMD) 
and a novel approach (4) Sparse Decomposition. The last two methods used 
are rather new and their potential use in time series analysis needs still to 
be explored. All times series were analysed in the same way with the four 
methods. After a short description of the methods we present the results and 
compare the results with respect to their capability of signal decomposition. 

1.2. Data 

We used four data sets to test the four methods for time series decompo- 
sition. The first data set (SYN) is a synthetic data set which is composed of 
a tidal signal calculated with a published Matlab®-script (Pawlowicz et al., 
2002), a long-term ramp function and two short-period events. In addi- 
tion noise with an RMS of 49.16 Pa is added. A step in the middle of 
the time series simulates a tectonic or oceanographic event. The second 
time series (MAR) consists of sea floor pressure observations from the Lo- 
gatchev Hydrothermal Field, located at the Mid- Atlantic Ridge (Gennerich 
and Villinger, 2011). Sea floor pressure measurements at this site were made 
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in order to monitor the magmatic and hydrothermal activity of the field which 
might express themselves as changes in sea floor pressure either due to sub- 
sidence or uplift or by tremor-like signals related to hydrothermal processes 
in the subsurface. The third and fourth records (GORKI and C0RK2) come 
from a so called CORK (Becker and Davis, 2005) installed at the sea floor 
off Vancouver Island. A CORK (Circulation Obviation Retroflt Kit) is a sea 
floor installation which is placed on top of a drill hole and isolates the hydro- 
logical - and pressure - regime inside the drill hole from the ocean. The main 
purpose of this installation is the long-term monitoring of pressure signals 
inside the drill hole in order to identify either slow earthquakes i.e. long pe- 
riod changes in the sub sea floor stress fleld or the response of the formation 
to seismic events which reveal the hydrological and permeability structure of 
the oceanic crust (e.g. Davis et al. (2011)). At every CORK installation sea 
floor pressures are recorded in addition to the downhole pressures to provide 
a local pressure reference signal. 

In all cases (MAR, CORKl and C0RK2) pressures were recorded with 
high-resolution absolute pressure sensors (http : / /www . paroscientif ic . com/ 
uwapp . htm) with a resolution of 7 Pa which is about 7 mm equivalent water 
hight change. More details on location and sample intervals can be found in 
table 1. 



Table 1: Details of the data sets used in the tests. 



Data set Start End Duration At Position 

[days] [min] Latitude Longitude 

SYN 41.7 60 48° 26'N 128° 43'W 

MAR 01.05.2008 31.05.2008 30 2 14° 45'N 44° 5' W 

CORKl 26.06.2003 15.09.2005 1046 10 



C0RK2 11.09.1996 15.09.2005 3292 60 



48° 26'N 128° 43'W 



2. Methods for Decomposing 

In this section we discuss the methods used to decompose the data. These 

are 

1. Harmonic Decomposition 



2. Wavelet Decomposition 

3. Empirical Mode Decomposition 

4. Sparse Decomposition 

The classical methods Harmonic and Wavelet Decomposition are well known. 
Therefore we have chosen them to have a comparison for the new methods 
Empirical Mode Decomposition and Sparse Decomposition. 

All methods have in common that they seek for a new reprensentation, 
as a superposition of several patterns, of the data, in most cases by basis 
transformation. 

2.1. Harmonic Decomposition 

The Harmonic Decomposition uses the Fourier Transform with its efficient 
implementation as Fast Fourier Transform (FFT). It is known since almost 
fifty years and used in many data processing routines. The basic idea is to 
transform the given time signal into the frequency space. In other words we 
are looking for a reprensentation with sinusiod pattern with different frequen- 
cies. We refer to Grochenig (2001) as a reference for a detailed introduction 
to Fourier Analysis. A much shorter introduction of the Fourier Transform 
is given in nearly every book about signal analysis, like Stark (2005). 

Decomposition of our time series with the help of Fourier Analysis starts 
with a FFT. As we know the frequencies and amplitudes of the tides from 
theoretical considerations and from observations, we can decompose the data 
into several components and remove the unwanted and deterministic ones. 
The remaining spectrum is transformed into a time series which contains 
only the non-tidal components. 

2.2. Wavelet Decomposition 

Another technique for decomposing the data is 
the Wavelet Decomposition which uses the Wavelet 
Transform. It is not as old as the Fourier Trans- 
form but it is also known very well. Like the Fourier 
Transform, it transfers the given time signal in an- 
other space. This time the basis vectors are scalings 
and translational displacements of a basic pattern, 
called mother wavelet, hence in this space the coef- 
ficients encode when an effect is measured and how ^^Sy® Wavelet dbAO 

(scaling function) 
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long this effect was. Tlierefore we can again decom- 
pose our data in features of a sliort and long period 

of time. Just like for the Fourier Transform there are fast implementations 
for the Wavelet Transform. For further informations about wavelets see Louis 
et al. (1997) or Stark (2005). 

First of all we have to choose the mother wavelet and a number of levels 
in which we want to decompose our signal. Unlike the Fourier Transform 
where the orthogonal functions are pre-defined, we have to choose the mother 
wavelet a priori. If we have extra knowledge about the data, we can improve 
the results by choosing a wavelet which is 'similar' to the data or rather the 
features we would like to extract. In our case the Wavelets Daubechies 10 - 
Daubechies 40 (often abbreviated by dblO and dbAO) fit the tidal components 
very well. After performing the Wavelet Transform with these wavelets, we 
obtained the best results by using dbAO, which is shown in figure 1. 

2.3. Empirical Mode Decomposition 

The third method we would like to introduce is the Empirical Mode De- 
composition and its enhancement the Ensemble Empirical Mode Decomposi- 
tion (EEMD). The EMD is a new method, invented by Huang et al. (1998), 
for decomposing a time signal. Since it is not as basic as the Fourier and 
Wavelet Transform, we present the algorithm here. The idea is to decompose 
the given time signal into a finite number of Intrinsic Mode Functions (IMF), 
which are functions which satisfy two conditions: 

1. The number of extrema and the number of zero crossings must differ 
at most by one. 

2. The mean value of the upper and lower envelope is zero at any time, 
where the upper envelope is computed by interpolating the maxima 
and the lower envelope by interpolating the minima. 

The task of the EMD is to find these IMF's and is done by the following 
sifting process, also shown in figure 2. 

1. Find all extrema of the given data. 

2. If the number of extrema is one or less, we have found all IMFs and 
terminate the algorithm. 

3. Compute the upper and lower envelope by interpolation the maxima 
and minima and the mean of these. 
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4. If the mean value is zero, we have found an IMF. Start the sifting 
process again with the data substracted by this IMF. Otherwise start 
the sifting process again with the data substracted by the mean. 



initialize 



substract mean 
no 



find extrcma 



(mcan^ 



yes 



T^extrema < 1 



no 



yes 



compute 
mean value 



substract IMF 



found IMF 



terminate 



Figure 2: The algorithm of the EMD. The coloured areas represent the loops of the 
algorithm. 

Also the EMD can be seen as a basis transformation, but this time the 
patterns, which are the IMFs, depend on the data. 

As we have seen we have done no a priori choices but choosing the in- 
terpolating scheme, which is predefined by cubic sphnes. Since abstracting 
an IMF reduces the number of extrema this algorithm terminates for every 
finite time signal. A disadvantage of this method is that it is empirical and 
has no solid theoretical foundation. Compared to the Fourier and Wavelet 
Transform the computing effort to decompose the data by EMD is higher 
but still tolerable. A detailed introduction to EMD is given in Huang et al. 
(1998). 

A challenge in decomposing time series is to prevent mode mixing. Phe- 
nomena of similar time scale should be in the same mode (here: IMF) and 
vice versa phenomena with a totally different scale are expected to have a 
different mode. For this purpose Wu and Huang suggested in 2009 a noise 
assisted data analysis method, called EEMD (Wu and Huang, 2009). For 
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a given signal white noise is added to the signal and then decomposed us- 
ing the EMD. At the end - after numerous trials - we take the mean of the 
IMF's. "By adding finite noise, the EEMD elimated largely the mode mix- 
ing problem and preserve physical uniqueness of decomposition. Therefore, 
the EEMD represents a major improvement of the EMD method." (Wu and 
Huang, 2009). Of course this is more time consuming than the usual EMD 
and we have to check if this is worthwhile in our application or not. 

2.4- Sparse Decomposition 

The last method we would like to present is the Sparse Decomposition. 
We assume that the given data is a superposition of some processes. Then 
we can find these by minimizing the residual 

||A;iXi + . . . + knXn - yh, 

where fci, . . . , fc„ are the possible patterns. Writing A;i, . . . , fc„ as the columns 
of a matrix K, then this is the same as 

\\Kx-y\\2. 

The Fourier and the Wavelet Transform can also be written as such a mini- 
mization problem by choosing the patterns as harmonics or wavelets respec- 
tivly. The advantage compared to the Fourier and Wavelet Transform is that 
we are not restricted to only one basis to find our decomposition, but we can 
incorporate several bases, which is called a dictionary. For our purposes 
we have chosen a dictionary with translational displacements and scalings 
of the pattern in figure 3, hence we combine the capability of the Fourier 
and Wavelet Transform. Furthermore, we can also add any other pattern to 
increase the physical meaning of the decomposition. 

We are looking for a solution which has a small residual and is sparse, 
i.e. it uses only a few patterns. A decomposition of the data into only a few 
patterns is likely to have physical meaning, especially when the patterns are 
chosen according to physical effects. This can be achieved by minimizing 

— y\\l + a\\x\\i, 

see Daubechies et al. (2004). The last term called penalty term, is the 

sum of the absolute values of x. The notion of Sparse Decomposition is also 
investigated in an applied manner in Chen et al. (1999). 
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u-^ ^ II — — ^ ^ ^ ii^^^ ^ ■ 

(a) harmonics (b) wavelets (c) peaks (d) blocks (e) linear trend 

Figure 3: The chosen dictionary which contains harmonic waves, a wavelet basis, a peak 
basis, blocks and a linear trend. 

For our application we have chosen the algorithm Regularized Feature 
Sign Search (RFSS), invented by Jin et al. (2009) and based on Lee et al. 
(2007), to minimize these functional. 

The choice of the dictionary is an advantage and disadvantage at the 
same time. If we have special knowledge of our data, like in most cases, we 
can improve the decomposition by choosing 'good' pattern. On the other 
hand an algorithm which requires many input parameters is often difficult to 
handle. A disadvantage is that the algorithm is much more demanding and 
complex compared to the other methods. We will see later on that this is 
worthwhile. 

Further information and the full algorithm is given by Jin et al. (2009). 
3. Results 

In this section we present and discuss our results. First of all we have to 
mention that the used data is preprocessed by subtracting the mean value, 
so the important features are emphasized. 

Our results for Harmonic and Wavelet Decomposition are obtained by us- 
ing Matlab®'s functions f ft and wavedec respectivly. For EMD and EEMD 
we have used the algorithms by Rilling (2007) and Wu (downloaded 04/2011). 
The RFSS implementation, used for the Sparse Decomposition, was accord- 
ing to Schiffler (2009). All algorithms can be found on the world wide web, 
see our references. 

The tidal extraction from the data sets is done very well by all algorithms, 
hence this does not distinguish the methods and is not discussed in the 
following. 




9 



3.1. Decomposition of SYN 

The results of the decomposition of the synthetic data are shown in figure 
4. Since it is a synthetic data set we can compare the computed decomposi- 
tion by the real one to evaluate their quality. 

As we see the methods' ability to decompose this data varies substantially. 
Harmonic and Wavelet Decomposition are able to decompose the data into 
features of three different types. The noise and the short effects are in the 
same component. The methods EMD and EEMD are not able to extract the 
high frequency features from the signal. All algorithms were able to separate 
the long-term effects from the rest but failed to distinguish between the trend 
and the step. 

Harmonic and Wavelet Decomposition as well as the Sparse Decomposi- 
tion correctly detect the short time features as well in time as in amplitude. 
Only large boundary effects in the results obtained by Harmonic and Wavelet 
degrade their quality. Methodically we can not obtain boundary effects in 
the results of the Sparse Decomposition, since we have chosen our dictionary 
to contain local and global pattern. Hence local features do not have to be 
reconstructed by global pattern. Additionally the RFSS is not based on filter- 
ing, which often causes boundary effects. Overall the Sparse Decomposition 
performs best in detecting the short time features. 

The last aim was to extract the long time effects. All algorithms perform 
quite well, but again the boundary effects degrade the Harmonic and Wavelet 
results. Sparse Decomposition is capable of extracting the trend and the step 
very well whereas EMD and EEMD reproduce only the trend satisfactorily. 

3.2. Decomposition of MAR 

The next data set to test is the one from the Mid-Atlantic Ridge and 
the results can be seen in figure 5. This time all algorithms provide a good 
decomposition into the short-, medium- and long-period features. Only the 
Sparse Decomposition is able to extract also the noise from the short-period 
components. There are just a few peaks visible in the short-period decom- 
position by Harmonic Decomposition and EMD which are probably aliased 
local eartquake signals. These are also extracted very well by the Sparse De- 
composition. The short-period features given by Wavelet show again some 
boundary effects. The decomposition by the EEMD has a basic noise level 
which is much too high. This is a result of the white noise added in every 
step of the algorithm but which is clearly not faded away with 100 runs of 
the EMD. 
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The trend and long-term feature extraction of the five methods are all 
nearly the same. In all cases pressures decrease initially over the first week 
and then increase by about 0.4 kPa over the rest of the time window. The 
Wavelet results differ from the others by significant boundary effects. 

3.3. Decomposition of GORKI 

The results of the third data set, GORKI, are given in figure 6. This time 
every method is again capable to decompose the data into the three main 
features. In this data set there are two short time effects whose amplitudes 
are a lot higher than the amplitudes of the other effects so we have also 
showed the short period effects zoomed in. 

The similarity between the short period features detected by Harmonic 
and Wavelet Decomposition is worth mentioning. The EMD and EEMD de- 
tect these features very similar in time but differ in detecting the amplitude. 
Also the EEMD is not able to detect short period features with a small am- 
plitude since the added noise prevents this. The Sparse Decomposition is 
able to detect these features as well and separates them additionally from 
the noise. 

The trends detected by all methods are again quite similar. Only the 
Wavelet trend has huge boundary effects. The behaviour of Harmonic De- 
composition at the boundary could be interpreted as a boundary effect, but 
we guess this is not the case, since Sparse Decomposition predict the same 
and is methodically unaffected by boundary effects. The same behaviour can 
be obtained from the EEMD's trend. 

3.4. Decomposition of C0RK2 

For the last data set the results are given in figure 7. Since this is a 
long data set and the sampling interval is only 60 min we focus on the long 
period effects. On the other hand we can not expect short time features 
with a shorter duration than 60 mm, like earthquakes, to appear in this data 
set, hence the short time components obtained from this data set has to be 
interpreted very carefully. 

Like with the dataset SYN, which sampling interval was also 60 mm, 
EMD and EEMD are not capable to decompose the features with higher 
frequencies than the tides. Contrary this is done very well by Harmonic, 
Wavelet and Sparse Decomposition. There are one huge and some smaller 
short period signals associated to earthquakes, which are again detracted 
from the noise by Sparse Decomposition. 
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The trends detected by all methods look quite similar. In any case we 
obtain first a downdrift of 3.5 kPa over 2 years, followed by an updrift the 
same amount in the next seven years. As seen at the other data sets the 
boundary effects are very high at the Wavelet Decomposition. 

3.5. Discussion 

A summary of the observations described in the previous sections is given 
in table 2. The evaluation of the different decomposition methods was based 
on the following more or less qualitative criteria: 

• How successful was the algorithm in decomposing the signal into a long 
and short period components? 

• Did the algorithm find a short period event in the presence of back- 
ground noise? 

• Does the algorithm create large effects at the boundaries of the time 
window? 

• How much CPU time is needed for the decomposition? 

From looking at table 2 one can conclude that none of the used methods 
for decomposition is perfect in every aspect. Harmonic and Wavelet are able 
to detect the short period features and perform also quite well in detracting 
the trend but their results, especially the results of Wavelet, show massive 
boundary effects. Because of this there is a high uncertainty if the shown 
decomposition is physically reasonable or not. On the other hand their com- 
puting efficiency is so high that the Harmonic and Wavelet Decomposition is 
nearly not restricted by the length of the given data. 

The performance of EMD and EEMD is directly opposite to Harmonic 
and Wavelet. The EMD and EEMD perform very well on detracting the 
tides and show only small influences of boundary effects but their capability 
of separating components with higher frequencies than the tides is restricted. 
At least when the data is recorded with a sampling interval of 60 minutes. 
Especially the EEMD fails because the added noise does not fade away in our 
experiments. Overall we have seen in every data set that the large additional 
computational effort of the EEMD compared to the EMD is not worthwhile 
for the analysis of sea floor pressure data. 

Very impressive is the capability of the Sparse Decomposition. It was able 
to detract the short period components as well as the large period components 
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Table 2: Overview of our results, indicates good, fair, and bad results. 
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in every case. Gontrary to the other methods, it was also able to detract the 
noise from the meaningful short period components. A hugh advantage is 
also that the Sparse Decomposition does not have any boundary effects. Only 
the computing efficiency needs to be improved. 
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Figure 4: Decomposition of data set SYN. There are shown from the top to the bottom 
the short, medium and long time components. The real decomposition can be seen in 
figure 4(f). EMD and EEMD were not capable to detect features with a higher frequency 
than the tides. 
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(a) Harmonic Decomposition 



(b) Wavelet Decomposition 
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(e) Sparse Decomposition 



Figure 5: Decomposition of the data set MAR. There are shown from the top to the 
bottom the short, medium and long time components. 
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(e) Sparse Decomposition 



Figure 6: Decomposition of the data set GORKI. There are shown from the top to the 
bottom the short, medium and long time components. 
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(e) Sparse Decomposition 



Figure 7: Decomposition of the data set C0RK2. There are shown from the top to the 
bottom the short, medium and long time components. Since the short time components 
has features of different amplitudes these are also showed zoomed in. EMD and EEMD 
failed to detect the short period features. 
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4. Conclusion 

To conclude our investigation we have to state that no method is perfect 
for every time series. For our data sets we got the best results by using the 
Sparse Decomposition. But at the moment it is not rccommcndablc for time 
series with more than 200,000 measurements. For this kind of time series 
we recommend the EMD for long period contributions and the Harmonic 
Decomposition for the short ones. EEMD fails to improve the EMD in de- 
tecting and is not recommended due to ist large computational needs. We 
obtained sometimes good results using the Wavelet Decomposition but they 
were never better than the results obtained by the Harmonic Decomposition, 
even in detecting earthquakes which is classically done by Wavelets. 

It is obvious that the comparison of the computing efficiency is somewhat 
unfair as algorithms to calculate the FFT and Wavelet Transform are highly 
optimized whereas algorithms for EMD, EEMD and especially RFSS are 
not yet written for speed. However it is very likely that FFT and Wavelet 
Transform will always be faster however the difference to EMD, EEMD and 
RFSS will most likely shrink in the near future due to the improvement of 
the algorithms and maybe by developing parallel algorithms for them. 
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